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tension  cut-off  capability  would  seem  to  most  easily  satisfy 
the  soil-structure  interaction  requirement.  However,  such  a 
formulation  does  not  provide  for  more  general  requirements 
involving  contact  conditions  of  structural  components.  Thus 
the  report  focuses  on  constraint-type  methods  in  addressing  the 
spatial  aspects  of  the  problem  of  implementing  dynamic  contact 
behavior  into  large-scale  finite  element  calculations  for 
civil -structural  systems. 

Of  the  constraint-type  methods,  one  in  particular  known 
as  the  augmented  Lagrangian  formulation,  is  deemed  most 
appropriate  to  the  specified  class  of  problems.  This  method 
is  a  combination  of  the  Lagrange  multiplier  and  penalty  form¬ 
ulations  and  aims  at  retaining  only  the  best  characteristics 
of  both  methods,  v^it  employs  an  augmentation  parameter  to 
vary  the  relati ve^participation  of  the  Lagrange  multipliers 
and  the  penalty  parameter  in  the  enforcement  of  contact 
conditions. 

Temporal  aspects  6f  the  problems  are  discussed.  Specifi¬ 
cally,  the  issue  is  raised  of  whether  or  not  special  handling 
techniques  for  temporal  integration  such  as  the  introduction 
of  impact- re  lease  conditions  are  required. 

Limited  numerical  experiments  suggest  that  very  little  addi¬ 
tional  computational  effort  is  required  to  obtain  the  full 
accuracy  of  the  Lagrange  multiplier  method.  Further,  the  work¬ 
able  range  of  the  penalty  parameter  is  shown  to  be  extended 
compared  with  the  conventional  method. 
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Methods  for  Solution  of  Dynamic  Contact  Problems 

by 

Robert  L.  Taylor,  PhD. 


1.  INTRODUCTION 

The  solution  of  contact  problems  by  the  finite  element 
method  has  received  considerable  attention  in  recent  years  [2- 
5,8,11-15,18,19,24-28,32-37,391.  In  general  the  problem  is  non¬ 
linear  and,  thus,  even  for  problems  composed  of  linear  elastic 
materials  is  non-trivial  to  solve.  For  transient  problems,  such 
as  those  of  interest  to  NCEL,  two  aspects  of  the  problem  need  to 
be  considered:  a.)  the  spatial  aspects  of  solution,  and  b.)  the 
temporal  aspects  of  solution.  For  the  spatial  part  of  the 
problem  solution  schemes  have  been  proposed  and  used  as  part  of 
finite  element  solution  algorithms.  Some  of  the  methods 
suggested  are  use  of  so  called  "gap"  elements  and  Lagrange 
multiplier  methods  [11],  penalty  function  methods  [7,12,13,24], 
perturbed  Lagrangian  methods  [36],  and  recently  augmented 
Lagrangian  methods  [37].  The  temporal  aspect  of  the  problem  has 
received  little  attention  in  the  literature.  In  the  work  of 
Hughes,  et.  al.,  [11]  it  was  stated  that  use  of  "impact-release" 
conditions  are  necessary  to  achieve  high  accuracy  in  the 
numerical  results.  This  is  contrasted  by  the  development  of 
Hallquist  in  which  no  such  conditions  are  included  yet  plausible 
results  are  achieved  [8].  Further  evaluation  is  necessary  before 


such  results  may  be  widely  applied. 

In  the  present  study  the  state-of-the-art  for  solving  large 


dynamic  contact  problems  modeled  by  the  finite  element  method  has 
been  reviewed.  The  problems  of  interest  to  NCEL  involve  such 
applications  as  large  dry-dock  facilities  which  may  experience 
strong  seismic  motions,  as  well  as  possible  accidental  blast 
loading  on  facilities.  Such  problems  are  characterized  by  the 
presence  of  materials  which  exhibit  very  strong  nonlinear 
responses  (e.g.,  soils,  concrete,  etc.).  The  anticipated 
deformations  for  many  of  the  conditions  anticipated  are  expected 
to  remain  small  so  that  linear  strain  measures  suffice.  The 
interaction  between  the  "structure”  (e.g.,  a  dry  dock  wall)  and 
the  surrounding  geological  media  and  the  interaction  of  system 
elements  (e.g.,  joints,  doors,  gates,  etc.)  are  expected  to 
involve  dynamic  contact  conditions  for  the  loadings  considered. 
Accordingly,  the  present  study  is  directed  to  possible  schemes 
for  including  dynamic  contact  modeling,  as  well  as  the  usual 
material  and  geometric  nonlinearity,  in  simulations. 

Past  studies  by  NCEL  have  focused  on  development  of 
appropriate  structural  modeling  capabilities  and  constitutive 
equations  for  the  systems  to  be  considered.  Constitutive 
equations  for  geologic  and  concrete  materials  have  been  developed 
and  studied  [9,10],  and  structural  modeling  methods  for 
plate/shell  systems  have  been  addressed  [16,31].  In  addition 
studies  have  been  conducted  by  Nour-Omid  and  Taylor 
[20,21,22,23,34]  on  effective  methods  to  solve  the  large 
algebraic  system  of  equations  resulting  in  this  class  of  non¬ 
linear  finite  element  problems.  Finally,  Shugar  and  Cox  have 
addressed  issues  related  to  combining  finite  element  and  boundary 
element  methods  [30].  In  future  efforts  NCEL  will  be  directing 
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efforts  to  bring  these  studies  together  in  a  design-oriented 
software  system  which  will  permit  engineers  to  make  quantitative 
assessments  of  the  behavior  of  the  soil-structural  systems 
subjected  to  transient  design  loads.  This  software  system  will 
provide  information  to  redesign  the  structure  for  loadings  which 
do  not  achieve  specified  limits  of  safety.  In  order  to  bring 
together  the  previous  efforts  it  is  necessary  at  this  time  to 
investigate  the  most  feasible  methods  for  including  the  contact 
interaction  with  the  other  salient  features  of  the  non-linear 
problems  involved. 


2.  METHODS  SURVEYED 

In  the  present  study  several  methods  were  surveyed  as 
possible  candidates  for  including  dynamic  contact  interaction 
effects.  The  methods  considered  include: 

a. )  Frequency  domain  methods, 

b. )  Use  of  non-linear  material  characteristics  and  no 

contact  models. 

c. )  Lagrange  multiplier,  penalty  and  other  associated 

methods , 

Each  of  these  methods  is  summarized  and  evaluated  below. 

2.1  Frequency  Domain  Methods 

In  the  analysis  of  seismically  loaded  soil-structure  systems 
a  frequency  domain  solution  is  often  advocated.  For  situations 
involving  non-linear  material  response  an  equivalent  linear 
material  model  based  upon  some  mean  measures  of  the  response  is 
introduced.  Such  methods  are  widely  used  in  the  design  of 
nuclear  power  plants.  Recently,  Roesset  has  mentioned  the 
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possibility  of  including  the  interaction  effects  of  separating 
boundaries  by  an  impedance  compensated  solution  [29],  No 
numerical  results  have  been  included  to  indicate  possible  results 
with  this  method.  The  primary  difficulty  with  use  of  frequency 
domain  solutions  is  the  development  of  effective  "mean" 
properties  for  each  of  the  model  constituents.  While  some 
success  has  been  achieved  for  analysis  of  soils  there  is  no 
indication  that  the  method  can  be  applied  in  general  to  a  wide 
class  of  problems.  Indeed  if  the  method  is  to  be  adopted  it  will 
be  necessary  to  qualify  the  results  by  an  extensive  set  of 
experiments  or  by  a  refined  objective  numerical  modeling  using  a 
more  accurate  deterministic  model.  Accordingly,  it  is 
recommended  that  this  approach  not  be  followed  at  this  time  and 
that  efforts  be  devoted  to  the  implementation  of  a  solution 
methodology  in  the  time  domain.  There  are  additional  reasons  to 
do  this.  For  strong  blast  loading  there  may  be  a  zone  in  which 
finite  deformation  effects  are  required.  Also  for  any  situation 
in  which  "unstable"  mechanisms  (e.g.,  buckling)  may  develop  it 
also  will  be  necessary  to  use  a  geometrically  nonlinear 
formulation^  At  this  time  there  is  no  evidence  that  frequency 
domain  solutions  can  incorporate  both  geometric  and  material  non¬ 
linearity  in  the  same  solution  procedure.  It  is  believed  that 
considerable  research  would  be  involved  to  even  assess  the 
feasibility  of  this  approach.  Finally,  the  frequency  domain 
approach  often  does  not  result  in  significant  savings  in  compute 
times.  On  the  other  hand  there  is  an  extensive  literature  on 
formulation  and  solution  of  geometrically  and  materially 
nonlinear  problems  in  the  time  domain. 


2.2  Use  of  Nonlinear  Material  Models  for  Interaction  Simulation. 

Several  non-linear  finite  element  computer  programs  have 

already  been  written  to  solve  specific  problems  which  involve 
both  geometric  and  material  nonlinearity.  The  most  notable 
programs  are  probably  ABACUS,  ADINA,  MSC/NASTRAN  and  the  DYNA  and 
NIKE  programs  developed  by  Hallquist  at  LLNL.  It  should  be  noted 
that  none  of  these  programs  is  a  design  oriented  system  such  as 
that  proposed  by  NCEL.  Many  of  these  nonlinear  finite  element 
programs  include  material  models  which  incorporate  a  "tension 
cut-off"  on  stresses  or  deformations.  In  modeling  granular 
materials  such  as  soils  the  inclusion  of  a  tension  limiting 
constitutive  assumption  is  essential  for  proper  modeling  of  the 
material  behavior.  In  modeling  soil-structure  interactions  the 
limiting  tension  will  effectively  provide  a  release  mechanism  of 
the  structure  from  the  soil  media.  In  this  class  of  problems 
there  may  be  little  need  to  incorporate  a  general  contact 
condition  in  the  modeling  of  the  system  since  the  tension 
limiting  aspect  of  the  soil  constitutive  equation  will 
essentially  eliminate  boundary  interaction  effects.  On  the 
otherhand,  for  interactions  between  structural  elements  there  is 
still  a  need  to  model  the  interaction  effects  as  reflected  by  a 
general  finite  element  contact  simulation. 

2.3  Lagrange  Multiplier  and  Penalty  Methods. 

A  methodology  to  handle  dynamic  contact  simulations  is 
included  in  both  DYNA  (an  explicit  method  finite  element  program) 
and  NIKE  (an  implicit  method  finite  element  program).  These 


programs  are  widely  used  at  LLNL  with  considerable  success  in 
modeling  dynamic  contact  interactions.  As  noted  previously, 
neither  of  these  programs  include  any  special  treatment  for 
"contact-release"  conditions. 

The  modeling  for  the  spatial  treatment  of  contact  problems 
may  be  incorporated  in  a  finite  element  program  by  introducing 
into  the  weak  form  of  the  balance  of  momentum  a  Lagrange 
multiplier  constraint  on  the  "gap".  By  monitoring  the  "gap" 
between  two  surfaces  the  Lagrange  multiplier  can  be  controlled  to 
introduce  the  correct  interaction  forces  whenever  a  contact 
condition  is  detected  (e.g.,  see  [11]).  In  a  discrete  form  of 
the  condition  it  is  necessary  to  describe  the  "gap"  in  some 
consistent  way.  Various  approaches  have  been  proposed  to  handle 
this  aspect.  In  the  earliest  works  deformations  were  considered 
as  small  and  it  was  assumed  that  nodes  on  one  surface  would 
"contact"  at  nodes  on  the  other  surface  (node-on-node  contact). 
This  leads  to  a  very  simple  form  of  the  Lagrange  multiplier 
method  to  describe  the  "gap".  Selecting  a  sign  convention  to 
which  a  positive  gap  corresponds  to  no  interaction,  a  negative 
gap  implies  penetration  and  it  is  merely  necessary  to  introduce 
the  Lagrange  multiplier  constraint  to  force  the  gap  to  a  "zero" 
value.  The  interaction  forces  are  "nodal"  in  nature  and  may  be 
monitored  to  ensure  that  tension  does  not  exist  on  the  surface. 
Accordingly,  the  form  of  the  Lagrange  multiplier  to  be  added  to 
the  weak  form  may  be  deduced  by  an  appropriate  linearization  of 
the  term 

(1)  -  F  (xx  -  x 2 ) 

where  F  is  the  nodal  interaction  force  (Lagrange  multiplier),  and 
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x^  and  X2  are  the  deformed  positions  of  a  "node"  for  surfaces  1 
and  2.  This  approach  is  extremely  simple  to  program  and  may  be 
included  as  an  "element"  subprogram  in  any  nonlinear  finite 
element  program  (e.g.,  see  Appendix  A  where  such  a  routine  is 
given  for  the  FEAP  program  developed  by  the  author  [31]).  While 
this  approach  is  very  simple  to  program  it  is  not  sufficiently 
general  to  permit  the  solution  of  realistic  contact  problems. 
Indeed  any  problem  involving  rolling,  sliding,  or  large 
deformations  will  invariably  lead  to  nodes  on  one  surface  not 
contacting  nodes  on  the  other  surface.  In  addition  modeling  of 
a  structure  for  dynamic  analysis  may  require  different  mesh 
spacing  for  each  component  or  different  material. 

It  is  relatively  straight  forward  to  extend  the  nodal 
Lagrange  multiplier  method  to  a  method  where  a  node  on  one 
surface  interacts  with  a  general  (non-nodal)  point  on  the  other 
surface.  For  two-dimensional  problems  a  Lagrange  multiplier 
term  may  be  introduced  in  the  form 
(2)  F  (  xx  -  (1-a)  x2i  -  a  x2j  ) 

where  "i"  and  "j"  define  two  nodes  on  the  2-surface  being 
contacted  by  the  node  on  the  1-surface  and  "a"  defines  the 
position  parameter  of  the  contact  point.  The  algorithm  is 
complicated  now  by  the  need  to  develop  search  methods  to  find  the 
nodal  pair  "i"  and  "j"  where  contact  occurs  and  to  determine  the 
value  of  "a". 

The  method  described  in  the  previous  paragraph  is  quite 
general;  however ,  some  deficiencies  have  been  noted  in  its 
application.  The  first  problem  is  that  there  is  a  bias  to  the 


treatment  of  the  surfaces.  Only  nodes  on  one  surface  are  in 
contact  with  the  other  surface.  The  nodes  of  the  other  surface 
may  not  be  in  contact  or  may  penetrate.  Indeed  the  second 
deficiency  is  this  latter  point.  If  several  nodes  exist  on 
surface  "2"  between  adjacent  contacting  nodes  on  surface  "1"  they 
may  penetrate  the  "contacting"  surface  in  an  uncontrolled  manner 
(see  Figure  1.).  This  is  most  undesirable  for  dynamic  problems 
where  inertial  effects  will  not  be  correctly  modeled  near  the 
contacting  surfaces.  Hallquist  proposed  to  treat  the  problem  by 
a  two  pass  approach.  In  the  second  pass  the  role  of  bodies  "1" 
and  "2"  are  interchanged  thus  averaging  the  contact  conditions. 
For  sufficiently  fine  meshing  this  is  an  effegtive  approach. 
However,  as  shown  by  Wriggers,  et.  al.,  for  coarse  mesh  regions 
the  results  may  not  be  esthetically  pleasing  [37].  Wriggers 
proposes  using  the  average  gap  on  segments  instead  of  nodal 
values  to  define  the  contact  conditions.  While  this  leads  to 
more  reasonable  pictorial  forms  of  the  contact  it  is  more 
difficult  to  determine  all  permutations  for  the  possible 
segments,  especially  for  3-dimensional  situations.  Some 
additional  effort  is  needed  to  define  appropriate  discrete  forms 
of  the  gap  for  3-dimensional  applications. 

An  additional  complication  from  the  introduction  of  a 
Lagrange  multiplier  treatment  of  the  contact  condition  is  that 
the  linearized  equations  of  the  finite  element  model  are 
indefinite.  Indeed  in  a  treatment  by  Newton  methods  with  a 
direct  solution  of  the  linear  algebraic  equations,  zero  diagonals 
occur  for  each  interaction  force  degree-of-f reedom.  For  static 
problems  direct  solution  procedures  without  pivots  can  fail.  It 
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has  been  traditional  to  not  use  pivots  when  solving  the  large 
linear  algebraic  equations,  hence  this  aspect  will  be  troublesome 
for  programs  using  direct  solution  procedures.  Often  the  problem 
may  be  avoided  by  a  judicious  numbering  of  the  degrees-of-freedom 
in  the  global  equations.  Since  many  programs  use  automatic 
renumbering  algorithms  it  will  be  necessary  to  introduce 
additional  checks  to  avoid  producing  a  zero  pivot  in  the 
renumbered  equations.  A  second  problem  is  that  additional 
degrees  of  freedom  are  introduced  into  the  global  equations  to 
account  for  the  interaction  effects.  Thus  even  if  direct 
solutions  are  replaced  by  appropriate  iterative  procedures  (e.g., 
the  Newton-Lanczos  method  recommended  and  studied  previously  for 
NCEL  [20,34])  these  additional  equations  must  be  retained. 
Finally,  checking  convergence  requires  separate  treatments  for 
the  force  and  displacement  parameters  to  properly  monitor 
performance.  It  is  possible  to  eliminate  the  interaction  forces 
"F"  from  the  global  equation  set  by  introducing  a  penalty  form 
for  the  Lagrange  multiplier  term.  This  is  very  effective,  and 
furthermore  removes  the  indefinite  character  of  the  resulting 
linearized  form  for  static  applications  of  stable  material 
models.  The  primary  difficulty  with  the  penalty  method  is  the 
need  to  select  an  appropriate  "penalty  parameter".  If  the 
parameter  is  too  small  undesirable  penetrations  will  result, 
whereas,  if  the  parameter  is  too  large  numerical  ill  conditioning 
problems  result.  Felippa  [7,8]  has  published  results  which 


illustrate  the  difficulties  which  may  result  from  using  a  penalty 
approach.  The  difficulties  he  cites  have  been  encountered  often 


by  us  in  using  penalty  methods  for  the  contact  problem  as  well  as 
in  the  enforcement  of  near  incompressibilty  conditions  in  finite 
deformation  elasticity  and  plasticity  problems  solved  by  the 
finite  element  method.  Precise  guidelines  cannot  as  yet  be 
established  to  indicate  when  the  selection  of  the  penalty 
parameter  will  lead  to  numerical  ill-conditioning  or  a  decrease 
in  rate  of  convergence  of  Newton's  method.  Accordingly,  it  is  a 
problem  which  must  be  addressed  further  before  large  finite 
element  programs  can  be  developed  for  general  application.  As  a 
pilot  study  an  effort  was  initiated  to  evaluate  the  effectiveness 
of  using  an  augmented  Lagrangian  approach  in  conjunction  with  a 
penalty  form  of  the  contact  constraint. 


3.  AUGMENTED  LAGRANGIAN  FORMULATIONS 

In  an  attempt  to  avoid  the  numerical  ill  conditioning 
problem,  the  Lagrange  multiplier  method  may  be  augmented  by  a 
penalty  term.  This  results  in  a  problem  formulation  which 
combines  the  best  features  of  both  methods.  The  global  equations 
are  normally  treated  using  the  penalty  term  in  the  linearized 
coefficient  matrix;  the  Lagrange  multiplier  term  is  included  as 
an  iterative  correction  to  the  right  hand  side  of  the  equations. 
The  method  is  generally  far  less  sensitive  to  the  penalty 
parameter  selected.  In  fact  the  penalty  parameter  may  be 
selected  as  a  small  number  and  iterative  improvements  are  used  to 
obtain  high  precision  satisfaction  of  the  constraint  term.  The 
number  of  iterative  corrections  required  to  achieve  convergence 
to  a  specified  tolerance  is  affected  by  the  size  of  the  penalty 
parameter  selected.  Thus,  the  algorithm  when  properly 
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implemented  requires  a  penalty  parameter  of  sufficient  size  to 
limit  the  iterations  to  a  reasonable  number.  The  proper  value 
for  the  penalty  parameter  in  an  augmented  Lagrangian  form  is  much 
smaller  than  that  required  for  the  penalty  term  alone.  For 
example,  a  problem  solved  by  the  penalty  method  and  requiring  a 
penalty  parameter  of  10^  to  achieve  "satisfaction"  of  the 
constraint  may  be  solved  using  a  value  of  102  to  102  for  the 
augmented  method  which  also  will  achieve  better  accuracy  with 
little  increased  effort. 

To  illustrate  a  use  of  the  augmented  Lagrangian  method  let 
us  consider  a  node-on-node  contact  formulation  as  defined  by  Eq. 
(1).  The  term  is  augmented  by  the  penalty  term  as  given  by  Eq. 

3. 

(3)  -  F  (  Xl  -  x2  )  +  0.5  (1  -  s)/E  (  F  )2 

The  penalty  parameter  is  given  by  "E"  and  "s"  is  a  parameter 
which  when  zero  defines  the  penalty  method  and  when  unity  gives 
the  Lagrange  multiplier  method.  Linearization  of  this  term  leads 
to  the  form 
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Felippa  has  discussed  several  possible  iterative  methods  to 
incorporate  the  Lagrange  multiplier  term  in  an  iterative  manner. 


Of  those  considered  the  hybrid  method  which  is  identical  to  an 
augmented  Lagrange  method  proposed  in  [37]  is  recommended.  Other 
methods  are  also  given  in  [1,17].  If  we  define  the  gap  by 

(6)  g(x)  =  xx  -  x2 

and  an  iterative  increment  to  the  position  by 

(7)  x(i+1)  =  x<i}  +  u(i) 

for  each  surface  then  the  Newton  equations  to  be  solved  at  each 
iteration  may  be  written  as 
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where  K^,  K2  are  the  tangent  stiffness  quantities  and  f^,  f2  are 


-E  I  I  Ut  | 

II  I 
(K2  +  E)  I  I  u2  I 


the  force  terms  for  bodies  1  and  2,  respectively;  and  "i"  is  the 
iteration  number.  As  suggested  in  [37]  the  force  F  may  be 
updated  using 

(9)  F(i+1)  =  F(i)  -  min  (  E  g(x)  ,  F(i)  ) 

As  the  method  converges,  g(x)  tends  to  zero  while  F  approaches 
the  solution  of  the  Lagrange  multiplier  method.  A  generalization 
of  this  algorithm  for  general  two-dimensional  contact  situations 
has  been  incorporated  into  the  finite  element  program  FEAP  using 
the  method  outlined  in  [37]  and  sample  problems  have  been  solved. 
Results  are  summarized  below. 


4.  NUMERICAL  EXAMPLES 

To  assess  the  performance  of  the  augmented  Lagrangian  form 
of  the  penalty  method  a  series  of  simple  test  problems  has  been 
analyzed.  The  only  method  assessed  under  the  current  effort  is 
the  perturbed  Lagrangian  form  of  the  problem.  In  our  implements- 


tion  it  is  not  possible  to  consider  an  arbitrary  meshing  of  the 
problem;  it  is  necessary  to  have  only  one  node  between  any  seg¬ 
ment  in  the  model.  This  is  not  a  limitation  of  the  method,  but 
merely  the  implementation  constructed  to  date.  As  will  be  seen 
in  the  recommendations  section  we  propose  to  use  the  two  pass 
penalty  method  for  the  three-dimensional  applications.  A  second 
reason  for  using  the  two-pass  algorithm  involves  performance  of 
the  perturbed  method  for  node-on-node  situations  (or  those  which 
are  node-on  node  to  a  numerical  tolerance)  where  the  algorithm 
has  some  aspects  of  ill-conditioning  in  defining  the  segment.  The 
perturbed  Lagrangian  does  have  an  advantage  as  will  be  illus¬ 
trated  in  one  of  the  examples  considered. 

4.1  Two  Spring  Constraint  Problem. 

The  first  example  considered  is  a  two  spring  model  for  which 
a  specified  gap  is  to  be  maintained  using  the  augmented  penalty 
method  described  above.  The  model  and  the  parameters  considered 
are  shown  in  Figure  2.  It  should  be  noted  that  the  "penalty 
parameter  will  be  placed  in  parallel  with  the  spring  "k2"  to 
maintain  the  gap  distance.  A  range  of  values  for  the  penalty 
parameter  was  considered  with  the  results  obtained  plotted  in 
Figure  3.  The  plateau  shown  in  Figure  4  denotes  maximum  number 
of  significant  digits  specified  for  defining  the  gap  and  is 
attained  in  very  few  iterations  (e.g.,  1-3  iterations  for  a  wide 
range  of  penalty  values).  This  example  illustrates  that  the 
augmented  method  requires  very  little  additional  effort  to  obtain 
the  full  accuracy  of  the  Lagrange  multiplier  method. 


4.2  "Rigid"  Punch  into  an  Elastic  Foundation. 

The  second  example  considers  a  stiff  punch  pressed  into  a 
flexible  elastic  foundation;  the  meshes  are  shown  in  Figure  4. 
This  example  has  been  considered  also  in  [37]  for  the  one-  and 
two-pass  penalty  algorithms  and  the  perturbed  Lagrangian 
algorithm.  In  this  example  the  modulus  of  the  punch  is  set  two 
orders  of  magnitude  larger  than  the  elastic  foundation  to 
represent  the  "rigid"  condition.  The  penalty  parameter  is  allowed 
to  vary  over  a  wide  range  of  values  ranging  from  10°  to  102®. 
The  accuracy  attained  (one  iteration  after  the  contact  region  is 
determined)  for  the  number  of  significant  digits  defining  the  gap 
is  shown  in  Figure  5.  The  maximum  accuracy  specified  is  8-digits 
and  is  obtained  for  the  augmented  method  for  all  values  of  the 
penalty  parameter  between  102  and  1016,  whereas  both  the 
perturbed  and  standard  (one  and  two-pass)  methods  attain  this 
accuracy  for  a  range  of  parameters  between  about  1010  and  10*4. 
The  converged  solution  accuracy  is  indicated  in  Figure  6.  This 
example  clearly  indicates  the  significant  advantages  of  the 
augmented  approach  over  a  conventional  penalty  method. 

4.3  Dynamic  Contact  of  an  Elastic  Block  on  a  Rigid  Surface 

The  next  example  considered  is  an  elastic  block  impacting  a 
rigid  surface.  This  is  a  dynamic  problem  which  has  a  wave  like 
solution  involving  a  situation  in  which  any  need  for  impact- 
release  conditions  should  be  evident.  The  mesh  and  properties 
are  shown  in  Figure  7.  The  deformed  mesh  (amplified  by  a  factor 
of  50)  is  shown  in  Figure  8  for  the  time  when  the  wave  traveled 
half  way  up  the  block.  The  block  had  a  release  time  which  was 
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near  the  theoretical  value.  The  stress  distribution  while 
qualitatively  correct  had  some  numerical  "noise”.  This  is 
probably  due  to  a  combination  of  factors  including  discretization 
error  of  the  mesh,  time  integration  errors  (an  implicit  Newmark 
implementation  was  used),  and  lack  of  impact-release  conditions. 
It  is  believed  that  the  latter  are  small  compared  to  the  other 
errors;  however,  additional  study  is  still  necessary  to  confirm 
this  impression. 

4.4  static  and  Dynamic  Contact  of  a  Sphere  on  a  Rigid  Surface 

The  last  example  considered  is  an  elastic  sphere  contacting 
a  rigid  surface.  The  mesh  and  properties  are  shown  in  Figure  9. 
The  results  for  three  static  loading  conditions  are  shown  in 
Figure  10.  For  the  mesh  considered  these  compare  very  favorably 
with  previously  published  finite  element  results  and  the  Hertz 
solution.  It  should  be  noted  that  the  perturbed  method  gives 
accurate  results  for  situations  where  only  three  segments  are  in 
contact.  The  dynamic  results  are  also  in  agreement  with 
previously  published  results  although  there  is  a  difference 
between  the  finite  element  results  and  those  of  Hertz  which  are 
based  upon  a  Ritz  solution  using  the  static  mode.  In  this 
example  the  load-time  curve  is  smooth,  appearing  similar  to  a 
sine-squared  loading,  and  the  impact-release  conditions  play  a 
very  minor  role  in  defining  even  a  theoretical  solution. 

5.  FINDINGS  AND  RECOMMENDATIONS 

The  present  study  has  considered  the  problem  of  dynamic 
contacts  modeled  by  numerical  methods.  In  particular  an 
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augmented,  perturbed  Lagrangian  formulation  has  been  considered 
and  is  recommended  as  the  best  available  method  for  inclusion  in 
a  design  package  to  be  developed  by  NCEL. 

The  augmented  Lagrange  multiplier  method  is  an  effective 
method  of  analysis  for  including  contact  capabilities  in  a  finite 
element  analysis  program.  It  permits  a  very  wide  range  of 
penalty  parameters  which  effectively  can  be  selected  to  avoid 
loss  in  numerical  accuracy  as  well  as  avoiding  ill-conditioning 
of  the  numerical  formulation.  Conventional  penalty  methods  lose 
accuracy  when  the  penalty  parameter  is  too  small  and  become  ill- 
conditioned  when  the  parameter  is  too  large.  As  illustrated  in 
even  simple  examples  conventional  penalty  methods  require  use  of 
a  very  carefully  selected  penalty  parameter  to  attain  viable 
solutions.  For  more  complex  problems  we  have  experienced 
difficulties  in  selecting  an  appropriate  value  of  the  penalty 
parameter  while  still  avoiding  the  above  cited  problems. 

The  above  findings  are  based  upon  limited  analyses; 
accordingly  the  following  recommendations  are  made: 

1.  A  three-dimensional  implementation  of  the  one-  and  two-  pass 
augmented  penalty  method  combined  with  a  gap  measured  by  the 
perturbed  approach  should  be  developed. 

2.  A  detailed  evaluation  of  the  performance  on  dynamic  problems 
involving  friction  and  impact-release  requirements  of  the 
type  needed  by  NCEL  should  be  made. 

3.  The  penalty  contact  method  should  be  evaluated  in 
conjunction  with  models  incorporating  the  non-linear 
constitutive  models  for  concrete  and  geological  materials. 
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General  design  goals  for  the  class  of  problems  to  be 
considered  should  be  constructed. 

5.  The  complete  problem  methodology  to  be  used  in  meeting  the 
design  goals  identified  in  recommendation  4  should  be  tested 
on  sample  problems. 

The  development  of  a  3-dimensional  implementation  is 
necessary  to  assess  possible  complexities  which  may  occur  in  the 
general  setting.  The  generalization  of  the  perturbed  Lagrangian 
method  to  3-dimensions  is  not  trivial  in  specifying  all  possible 
segment  configurations.  Recently,  it  has  become  apparent  that 
measuring  the  average  gap  for  the  two-pass  method  is  nearly 
equivalent  to  the  perturbed  method.  Accordingly,  it  is 
recommended  to  consider  the  development  of  the  two-pass  method, 
but  with  an  average  gap  used,  instead  of  individual  gaps  at 
nodes,  to  define  the  contact  surface.  An  important  aspect  of 
using  Newton's  method  is  the  proper  derivation  of  a  tangent 
matrix.  Thus,  efforts  must  be  made  to  derive  the  consistent 
tangent  array  for  all  non-linear  terms  in  the  structural  model  - 
including  the  contact  terms. 

The  behavior  of  dynamic  contact  problems  analyzed  without 
use  of  special  contact-release  conditions  is  promising.  It  has 
become  apparent  in  our  recent  efforts  in  solving  a  variety  of 
non-linear  problems  that  care  must  be  exercised  in  developing  the 
algorithm  for  the  step-by-step  integration  of  the  equations  of 
motion.  A  methodology  is  now  used  in  FEAP  which  has  performed 
well  for  analyses  which  include  large  motions  (e.g.,  analyses  of 
free  flight  of  structures  representing  flexible  space-craft),  as 
well  as  those  which  have  non-linear  material  behavior.  Our 
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analysis  of  the  dynamic  contact  problem  has  also  used  this 
formulation.  To  date,  we  have  used  Newmark’s  method  with  the 
start  condition  at  each  time  step  the  last  converged 
"displacement"  state.  Initial  velocities  and  accelerations  are 
deduced  from  the  Newmark  equations  to  be  consistent  with  this 
starting  condition.  We  have  also  been  extremely  careful  to 
include  a  proper  linearization  of  specified  boundary  motions  so 
that  they  are  consistent  with  the  algorithm  for  the  active 
degrees  of  freedom.  The  implementation  is  then  completed  by 
computing  iterative  corrections  to  the  nodal  displacements  and, 
subsequently,  the  nodal  velocities  and  accelerations.  Non-linear 
material  models  are  integrated  at  each  "stress  point"  using  a 
local  Newton  iteration  for  the  full  increment  of  motion  computed 
in  the  time  step.  We  know  that  the  implementation  details  for 
the  Newmark  method  are  sensitive  in  solving  the  contact  problem. 
In  using  a  previous  implementation  it  was  necessary  to  introduce 
contact-release  conditions  in  order  to  obtain  viable  results  [2], 
Further  analyses  are  required  to  ensure  that  accurate  results  may 
be  obtained  using  the  current  implementation  on  a  wide  range  of 
problems.  In  addition  effort  should  be  devoted  to  developing  a 
formulation  which  includes  a  consistent  treatment  of  friction  in 
the  formulation.  Recent  applications  of  operator  split  methods 
have  been  used  to  clearly  explain  the  "radial  return"  algorithm 
for  solving  elasto-plasticity  problems.  The  similarity  between 
friction  problems  and  e 1  a s t i c -p 1  a s t i c  problems  cannot  be 
overlooked.  It  is  recommended  that  an  operator  split  formulation 
be  pursued  in  an  attempt  to  develop  a  consistent  tangent 
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formulation  for  frictional  contact  problems.  This  formulation 


will  need  to  be  evaluated  for  dynamic  problems  to  ensure  that 
impact-release  conditions  are  still  not  required. 

The  general  design  goals  and  loading  states  to  be  considered 
to  achieve  these  goals  needs  to  be  clearly  stated  by  NCEL  to 
ensure  that  the  additional  developments  are  not  required  in 
developing  a  design  package.  This  will  involve  some  preliminary 
analyses  to  assess  the  levels  of  deformation  and  stress  which  can 
be  anticipated  in  the  facilities  of  interest  to  NCEL.  It  is 
particularly  important  to  assess  whether  large  strain  effects  or 
instability  may  result  in  any  of  the  problems.  This  phase  will 
identify  the  salient  feature  which  the  analysis  program  must 
include  in  order  to  develop  a  design  oriented  program  system. 
Some  attention  needs  to  be  devoted  to  the  methods  of  showing 
results  to  the  designer.  For  example,  use  of  sensitivity 
analysis  may  provide  important  insight  for  the  designer  to  effect 
changes  in  the  system  parameters  and  thus  achieve  the  design 
objectives . 

The  above  recommendations  involve  considerable  effort  and 
must  be  be  given  a  priority  ranking.  For  example,  in 
establishing  the  design  goals  it  will  be  necessary  to  specify 
computing  requirements.  This  in  turn  depends  upon  the 
capabilities  of  particular  computers  available  to  NCEL  and  its 
contractors.  This  is  currently  a  rapidly  changing  area  and 
decisions  must  not  be  made  too  early  in  the  setting  of  goals.  On 
the  otherhand,  the  assessment  of  the  behavior  of  the  augmented 
Lagrangian  method  with  the  other  material  models  may  proceed 
independent  of  which  computer  will  ultimately  be  used. 


Accordingly,  while  effort  should  be  given  early  to  all  aspects  of 
the  problem  those  involving  additional  developments  should  be 
addressed  first. 
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Figure.  5.  Results  after  1  iteration  for  the  3  block  problem 
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APPENDIX  A.  Contact  Element  for  FEAP  -  Node-on-Node  Case 


subroutine  elmtl2<d,ul ,xl , i  x  ,  1 1 ,s,p,ndf,ndm,nst, i  sw> 
implicit  double  precision  <a-h,o-z> 
c....  contact  element  -for  feap 

dimension  d<l),xl<ndm,l),u!(ndf,l>,ix<l),tl<l),s<nst,l),p<nst) 
dimension  wd<  2) 

c ommon/c  da  ta/o , he ad  < 20 ) , n  umn p , n ume 1  , n  umma t , n e n , neq  , ipr 
common/e  1  data./ dm  , n  ,ma  ,mc  t ,  i  e  1  ,  ne  1 
save  /cdata/ ,/e  1  data/ 
data  wd/4hsl i p , 4h 1 ock/ 
c....  transfer  to  correct  processor 

go  to  <1, 2, 3, 3, 5, 3, 5, 2, 2, 2),  isw 
c....  input  the  material  properties 
1  re ad <5, 1000)  i df , i d i r , i st , d< 4) , i c e 1 m , d( ?> 

is  t  =  max  0 ( 1  , m i n  0 ( 2 , i s  t  >  > 
i f <d<4> . 1 e .0 .0d0)  d<  4)  =  l.e-5 

uri t e ( 6 , 2000  >  i  df , i d i r ,wd( i s t ) , d<  4 ) , i c e 1 m , d <  9) 

d<  8)  =  icelm 

d< 1 )  =  i  df 

d  <  2  >  =  i  d  i  r 

d  (.  3 )  =  i  s  t 

return 

....  check  element  for  errors 
continue 
re  turn 

....  compute  the  element  tangent  array 
i  df  =  d <  1  > 
idir=  d<  2) 
i  st  =  d(  3) 

qap  =  <  u  1  <  i  df  ,  3  >  -u  1  (  i  df  ,  1  )  >  +  ( x  1  (  i  d  i  r  ,  3 >  -x  1  (  i  d  i  r  ,  1  >  ) 
tt  =  ul < i df ,2) 
e  =  0 .0 

i f  < gap . 1 t . d  < 4 ) )  e  =  1.0 
i f <  tt . 1 t .-d<4) >  e  =  0.0 
i f  <  gap . 1 t . -d  <  4> )  e  -  1.0 
i f  < i sw . ne . 3)  go  to  4 
go  to  <30,31),  i s t 
....  sliding  contact 
0  i  =*  ndf  +  i  df 
s< i , i df >  =  e 
s<  i  df  ,  i  )  =  e 
s< i , i +ndf )  =  -e 
s< i +ndf , i )  =  -e 
s(i,i>  =  1.0  -  e 
r e  tur n 

c....  perfect  stick  contact 

31  do  32  i  =  1 , ndf 
s< i  , ndf + i >  =  e 
s<  ndf + i  , i )  =  e 

s< ndf + i , ndf + i )  =  1.  -  e 
s< ndf +ndf + i , ndf + i )  =  -e 

32  s< ndf ♦ i ,ndf +ndf + i )  =  -e 


return 

compute  and  output  the  element  variables 
continue 

i f ( i sw . eq . 6)  go  to  6 
wr i te < 6 , 2002)  n,ma,tt,e 
return 

compute  element  mass  arrays 

continue 

return 

compute  the  element  residual  vector 

p  <  i  df +ndf )  =  gap 

go  to  (60,61),  i st 

P ( i df  )  =  -e*t  t 

p ( i df +ndf +ndf >  =  e*tt 

if(e.eq.O.OdO)  p<idf+ndf>  =  -tt 

go  to  64 

do  62  i  =  1  ,  ndf 

p ( i )  =  -e*u 1 < i , 2) 

P  (  i  +  ndf  +  ndf  >  =  e*ul(i,2) 

•  -  (  e  .  i q  .  0  . 0 d0  :  n d-f  +  i  )  =  -ul  <  i  ,2) 
i  ce  1  m  =  d(  8) 

i  f  (  i  c  e  1  m .  eq  .  n  >  p<  i  of  ♦nd-f +  n  dr  =  o  :  of  ♦  nd-f  +  nd-f )  +  d<  ?> 
return 
f  orma  t  s 

f  orma  t(3i5,fl0.0,  i  5 ,  f  1 0 . 0  > 

f  orma  t  ( 5x  , '  con  tac  t  e  1  emen  t  '/l  Ox  , '  con  tac  t  d.o.-f.  '  ,i5/10x 
'coordinate  direct.' , i 5/ 

'contact  condition  '  ,a4/'10x  , '  tol  e  ranee  ',el5.5/ 
l  'forced  e 1 emen t ', i 6/" force  value  ',el3.5/> 

formats' el mt ' , i 5 , 'mat  1 ' , i  5,  ' force ' , e 1 2. 3 , ' e' ,  f  6 . 2> 


